#include "vtkMeshExport.hpp"
#include "vtkMeshImport.hpp"
#include "vtkFieldExport.hpp"
#include "vtkFieldImport.hpp"
#include "sweSolver.hpp"

#include "meshTypeFunctions.hpp"

#include <catch2/catch_test_macros.hpp>

#include <cmath>

using namespace hms;
using namespace hms::vtk;
namespace xml = tinyxml2;

TEST_CASE( "Domain field import" ){

	auto test = []( auto&& meshDummy){
		using MeshType = remove_qualifiers<decltype(meshDummy)>;

		MeshVariant meshVar { uniMeshAs<MeshType>(3,4,1.3,0.8) };

		const MeshType& mesh { std::get<MeshType>(meshVar) };

		Time time {10};
		time.writeTimeInterval(1);

		SWESolver<MeshType> solver {mesh, time};

		std::unique_ptr<xml::XMLDocument> doc { meshToVtk(meshVar) };
		xml::XMLElement* vtk { getVtkNode( doc.get() ) };
		writeDomainFields(vtk, solver);

		/* get CellData node, which holds the internal field values */
		xml::XMLElement* cellData { getChildNodeByName(
			getChildNodeByName(
				getMeshNode<MeshType>( getVtkNode( doc.get() ) ),
				"Piece"
			),
			"CellData"
		) };

		/* now we can write to individual internal fields */
		xml::XMLElement* vDataArray { findFieldData(cellData, "Velocity" ) };
		xml::XMLElement* dDataArray { findFieldData(cellData, "WaterDepth") };
		xml::XMLElement* zDataArray { findFieldData(cellData, "BottomElevation") };

		vDataArray->SetAttribute("isUniform", false);
		dDataArray->SetAttribute("isUniform", false);
		zDataArray->SetAttribute("isUniform", false);

		vDataArray->SetAttribute("isFunction", false);
		dDataArray->SetAttribute("isFunction", false);
		zDataArray->SetAttribute("isFunction", false);

		SECTION( "read non-uniform value lists" ){
			ArrayXXs vRand, dRand, zRand;
			vRand.resizeLike( solver.v().domainValues() );
			dRand.resizeLike( solver.d().domainValues() );
			zRand.resizeLike( solver.z().domainValues() );
			vRand.setRandom();
			dRand.setRandom();
			zRand.setRandom();

			matrixToText( vDataArray , vRand);
			matrixToText( dDataArray , dRand);
			matrixToText( zDataArray , zRand);

			readDomainFields(vtk, solver, false);

			REQUIRE( solver.v().domainValues().isApprox(vRand) );
			REQUIRE( solver.d().domainValues().isApprox(dRand) );
			REQUIRE( solver.z().domainValues().isApprox(zRand) );
		}
		SECTION( "read uniform values" ){
			Array2s v { 0.3, -1.4 };
			ArrayNs<1> d { 3.4 }, z { -9.9 };

			vDataArray->SetAttribute("isUniform", true);
			dDataArray->SetAttribute("isUniform", true);
			zDataArray->SetAttribute("isUniform", true);

			matrixToText( vDataArray , v);
			matrixToText( dDataArray , d);
			matrixToText( zDataArray , z);

			readDomainFields(vtk, solver, false);

			REQUIRE( solver.v().domainValues().isApprox(v.replicate(1,12)) );
			REQUIRE( solver.d().domainValues().isApprox(d.replicate(1,12)) );
			REQUIRE( solver.z().domainValues().isApprox(z.replicate(1,12)) );

			SECTION( "Skipping z field" ){
				solver.z().domainValues() = 16;
				readDomainFields(vtk, solver, true);
				REQUIRE( solver.z().domainValues().isApproxToConstant(16) );
			}
		}
		SECTION( "read non-uniform function(s)"){
			vDataArray->SetAttribute("isFunction", true);
			dDataArray->SetAttribute("isFunction", true);
			zDataArray->SetAttribute("isFunction", true);

			writeFunction(vDataArray, "x", 0);
			writeFunction(vDataArray, "y", 1);

			writeFunction(dDataArray, "z", 0);
			writeFunction(zDataArray, "sin(x*y)", 0);

			readDomainFields(vtk, solver, false);

			Index nDomainCells { mesh.nDomainCells() };

			Array2Xs vCompare (2, nDomainCells );
			ArrayXs dCompare ( nDomainCells );
			if constexpr ( std::is_base_of_v<StructMeshBase<MeshType>, MeshType> ){
				Index
					nx { mesh.size().x() },
					ny { mesh.size().y() };
				for ( Index j {0}; j<ny; ++j ){
					for ( Index i {0}; i<nx; ++i ){
						Index cell { j*nx+i };
						const auto& centroid { mesh.cell(i,j).centroid() };
						vCompare.col( cell ) = centroid;
						dCompare(cell) = std::sin( centroid.prod() );
					}
				}
			} else {
				for ( Index i{0}; i<nDomainCells; ++i ){
					const auto& centroid { mesh.cell(i).centroid() };
					vCompare.col(i) = centroid;
					dCompare(i) = std::sin( centroid.prod() );
				}
			}

			REQUIRE( solver.v().domainValues().isApprox(vCompare) );
			REQUIRE( solver.d().domainValues().isApprox( dCompare.transpose() ) );
			REQUIRE( solver.z().domainValues().isApprox( dCompare.transpose() ) );
		}
		SECTION( "Read uniform functions" ){
			vDataArray->SetAttribute("isUniform", true);
			dDataArray->SetAttribute("isUniform", true);
			zDataArray->SetAttribute("isUniform", true);

			vDataArray->SetAttribute("isFunction", true);
			dDataArray->SetAttribute("isFunction", true);
			zDataArray->SetAttribute("isFunction", true);

			writeFunction(vDataArray, "sin(3/4*pi)", 0);
			writeFunction(vDataArray, "g/2", 1);

			writeFunction(dDataArray, "123", 0);
			writeFunction(zDataArray, "e", 0);

			std::cerr.setstate(std::ios_base::failbit);
			readDomainFields(vtk, solver, false);
			std::cerr.clear();

			Index nDomainCells { mesh.nDomainCells() };

			Array2s vCompare { std::sin(3./4.*math::pi), phys::g/2 };
			scalar dCompare { 123 }, zCompare { math::e };

			REQUIRE( solver.v().domainValues().isApprox(
				vCompare.replicate(1, nDomainCells)
			) );
			REQUIRE( solver.d().domainValues().isApproxToConstant( dCompare ) );
			REQUIRE( solver.z().domainValues().isApproxToConstant( zCompare ) );
		}

	};

	SECTION("Mesh"){
		test( Mesh{} );
	}
	SECTION("UniMesh"){
		test( UniMesh{} );
	}
	SECTION("RectMesh"){
		test( RectMesh{} );
	}
	SECTION("StructMesh"){
		test( StructMesh{} );
	}
}